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(N , ABSTRACT 

Q I This work assembles some basic theoretical elements on thermal equilibrium, stability 

conditions, and fluctuation theory in self- gravitating systems illustrated with a few examples. 
Thermodynamics deals with states that have settled down after sufficient time has gone by. 
CN I Time dependent phenomena are beyond the scope of this paper. While thermodynamics 

• is firmly rooted in statistical physics, equilibrium configurations, stability criteria and the 

T-H ■ destabilizing effect of fluctuations are all expressed in terms of thermodynamic functions. The 

^ I work is not a review paper but a pedagogical introduction which may interest theoreticians 

in astronomy and astrophysicists. It contains sufficient mathematical details for the reader 

■ to redo all calculations. References are only to seminal works or readable review sQ. Dehcate 
^ [ mathematical problems are mentioned but are not discussed in detail. 

(N : 

■ Keywords: Gravity - Thermodynamics of non-extensive systems - Statistical mechanics 
r~| . with long range forces. 
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■ 1 Introduction 
> . 

^ , A strident aspect that puts self-gravitating systems in a class of their own is that they can- 

^ ' not be small sub-systems of a large ensemble except in very imaginary situations. A small 

sub-system of a large self-gravitating system has its energy almost entirely determined by the 
large system. Moreover, the sum of energies of all the small sub-systems is not equal to the 
total energy of the system in contrast with a basic tenet of classical thermodynamics. This 
has far reaching consequences. Much of the rich harvest of classical and quantum thermody- 
namics as exposed in Landau and Lifshitz's Statistical Physics [^] is useless in gravitational 
thermodynamics. For instance stable and unstable isolated self-gravitating systems may have 
heat capacities of both signs. This is also the case for a system in a heat bath. Thus Cy > 
is no more a signature of stable thermal states like in classical thermodynamic ensembles. In 
this respect gravitational thermodynamics appears slightly exotic though not to astronomers 
as pointed out by Lynden-Bell |49| in his review on progress in understanding of gravita- 



tional thermodynamics. It is also what makes thermodynamics of self-gravitating systems 
interesting. 

The paper is based on the Gibbs density of states function in the meanfield approximation. 
Astronomers and astrophysicists are more familiar with Boltzmann entropy. Gibbs definition 



^Except in the ultimate sub-section. 
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is more general in that it provides Boltzmann entropy as a first approximation as well as 
stability conditions and fluctuations. It has a wider range of applications and extends to 
meanfield theory. The Gibbs approach is in our view the least frustrating way to understand 
statistical thermodynamics and provides the best logically connected approach to the subject. 

Section 1 begins with a review of statistical equilibrium theory in the simple case of 
N identical point masses in an external field first and gravitationally interacting particles 
in a finite volume next. The case = 2 is interesting and has been studied in detail 
by Padmanabhan |5^. Here we consider N >> 2 and situations in which particle pair 
formation has a negligible role on relevant time scales. In this case meanfield theory is a very 
good approximation. A steepest descent calculation gives equilibrium configurations of stable 
as well as metastable states. Metastable states are particularly interesting in gravitational 
thermodynamics^. The method provides also stability conditions. 

In section 2 we deal with thermodynamic functions and stability conditions in general. 
The latter are first examined with the help of Gibbs' density of states function but are 
subsequently translated into conditions involving only thermodynamic functions. The theory 
applies to isolated non-extensive systems, the type of system considered in section 1, but is 
of much broader applicability. It predicts in particular that under very general assumptions 
stable thermal equilibrium configurations of isolated systems near instability have always 
negative heat capacities which turn to be positive when systems become unstable! In classical 
thermodynamics a small thermally stable subsystem of an ensemble has always a positive heat 
capacity. These results imply that stable isolated systems near instability are necessarily 
unstable if put in a heat bath. This well known "nonequivalence of ensembles" has no 
analogue in classical thermodynamics. 

Section 3 develops the theory of fiuctuations. We deal with fiuctuations of temperature 
as a generic example. Fluctuations of other thermodynamic functions can be dealt with in 
a similar way. The theory has also broad applicability; again results are expressed in terms 
of thermodynamic functions and hold only near instability. For nonequivalent ensembles the 
behavior of equilibrium near instability departs considerably from classical systems and the 
effect of fiuctuations is consequently quite different. 

Section 4 gives two examples of application with some detail: isothermal spheres in various 
ensembles and liquid ellipsoids. The purpose here is to illustrate the power and limitations 
of thermodynamics in self- gravitating systems. The last sub-section reviews briefly various 
other applications in astronomy where thermodynamics has been used or might have been 
with considerable benefit. 

Each section has its own summary, conclusion and/or comments. 

2 Statistical thermodynamics 

It may interest the reader to see the explicit connection between the entropies of Gibbs and 
Boltzmann by starting with the simpler case of a system in an external field first. Here 
some of the mathematical techniques appear in a more transparent form. We go then over 
to self-gravitating systems in the meanfield approximation which is the main object of this 
section. 

(i) Gibbs's density of states function and Boltzmann entropy. 

^Most globular clusters j50| and rich clusters of galaxies have their cores in a quasi thermalized 
metastable state. 
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Consider first the simple case of identical point particles of mass m in a volume V, 
with perfectly reflecting massless walls in an external gravitational potential U{r). The total 
energy H, the sum of kinetic Hx and potential energy Hp of the system is a constant of 
motion. Let E be its value. The total energy H is a function of the coordinates of each 
particle in phase space (rj,pj), (i = 1,2,..., A^): 

H = Hk + Hp = Y, ^pi + ^ mllin) =Y,E, = E. (2.1) 

It I 

The total energy H is here equal to the sum of the energies of individual particles Ei. The 
density of states function N, m) or in short ^{E) is the sum of all possible states^ of 

the system with energy E divided by iV!: 

Q{E) = 6{E - H)Y[duJi where duji = (fricfpi (-00 < pi < +00); (2.2) 

i 

6{E — H) n diOi is known as the microcanonical distribution p^ ] and the basic justification 

i 

for this distribution is the ergodic hypothesis. The Gibbs entropy Sg is the logarithm of Vt{E) 
divided by Boltzmann's constant k. Here we take k = 1 which amounts mainly to measure 
the temperature in units of energy. Thus 

n{E)=e^^. (2.3) 

To evaluate the Gibbs entropy we first replace 5{E — H) by a Bromwich integral of 1 (its 
inverse Laplace transform) in brackets: 



n{E) = 



1 /-fe+ioo (3{E-Y,Ei) 

e - dp 



2Tri Jb 



-IQO 



\{duji. (2.4) 



This expression for J7 can be rewritten after a slight rearrangement of terms as 



L % J 

^, the quantity between brackets, is easily seen to reduce to 

^ = A^J^J e'^^^^du^j where ^0 = + rnU{r) and duj = d^rd^p. (2.6) 

Eq is the energy of of one particle with coordinates (r, p). We can of course also write ^ as 
follows 

N'=co -I / p _ \N' 

^= E ^ /^"^''"^^ (2-7) 

N'=0 • V"' / 



■^To avoid any confusion we use here a tilde over all quantities with the same symbols and the same meanings 
as in self-gravitating systems. 

^The Gibbs density fo states function has dimensions [Q] = erg^^"^ sec'^^ . Thus f2 like later Q is defined 
up to a constant that depends on the units chosen. 
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We then approximate (Jjv'at by 5{N' — N) since N >> 1 and replace immediately 6{N' — N) 
by a Bromwich integral like we did with 6{E — H): 



1 

2m J a 



a+ioo N'=oo 

y — 

^ ATM 

Af'=0 



N' 



(2i 



^ can now be rewritten in this form 

2 /-a+ioo 



2'7rz J 0—200 



/ = e' 



a-l3Eo 



and if we substitute this expression for ^ back into ( |2.5D we obtain 



1 



a+ioo /-b+ioo 



-ioo Jb—ioo 



a = -aN + pE+ / fduj. 



(2S 



(2.10) 



We shall now evaluate Q(E) by a steepest descent technique Q as follows. We first look for 
the extremum of a{a,f3); this is obtained from the following equalities: 



da 



da 



da 



— = -N + fduj = Q , — = E- Eofduj = 0. 



5/3 



(2.11) 



The equations define a point (de, Pe) in (Oj /?) space in terms of E and N and a corresponding 
equilibrium value fe for /: 



(2.12) 



/e is the Boltzmann distribution of energies in the external field U{r) calculated at the point 
of extremum. The extremal value ae of a is the Boltzmann entropy S and it is easily seen 
that 



S = ae. 



feHfe)dLJ + N. 



(2.13) 



We may thus write ( p. 10 ) as follows introducing 6^ a which represents the sum of all the terms 
of order higher than one in a Taylor expansion of a near the extremum ae- 



(27ri) 



a+ioo /"fe+ioo 



a— ioo Jb—ioo 



~^dadf3. 



(2.14) 



Since the exponent 5'^a is of order the integrant in ( 2.14| ) is very steep and a good approx- 
imation is obtained by limiting the integrations to terms of order two. The terms of order 
two are readily found; with 5a = a — ote and 5(3 = (3 — (3e we may write 



5^a 



1 
2 
1 
2 

2 



a 



da^ 



\dadl3 j 



{5af + 2 

N{5af - 2E5a5(3 + NE^{5(3f 
{5a - mpf + (^0 - E'of5p'^ 



5a5(3 + 



'd^ 



a 



dp^ 



{5(5f 



(2.15) 



in this formula mean values calculated in the phase space of a particle are denote with an 
overbar; for instance: 

1 



En 



N 



Eofeduj 



(2.16) 
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We can now fix the position of the imaginary axis for the complex variables (a, /3) by intro- 
ducing a new set of variables which vary between ±00: 



a* = i{5a — Eo6f3) and (3* = i6f3. 
This gives to ( |2.10| ) the following form: 



n{E) = e' 



e 2 



The integrals are easily found and the end product is this: 



1 



SgC^S- In(iV) - ln(27r) - - In 



(^0 - Eo)' 



(2.17) 



(2.18) 



(2.19) 



This shows incidentally that S is a maximum of a for any U; in other words the system is 
stable and ^1{E) is convergent. Moreover, since S is of order and >> 1 the last three 
terms are clearly negligible compared to S so that 



iSg — *S* . 



(2.20) 



This is the main result of this first subsection. 



(ii) The Gibbs density of states function for self-gravitating systems. 

If the same system of particles of mass m in a volume V are interacting gravitationally, 
the total conserved energy is 



2 " iFj — r,-| 



E. 



(2.21) 



G is the gravitational constant. The density of states function Q{E) is the same as ( |2.2| 
without tildes: 

n{E) = ^1 6iE - H)l[dcoi. (2.22 



As is well known, Q diverges. This can easily be seen as follows^. Take equation (|2.22D , 
replace again 5{E — H) by a Bromwich integral of 1 with variable E — H = E — {Hk + Hp)] 
in this way we obtain an expression similar to (|27 



rb+ioo 



ME-Hk-Hp) 



dp 



Y[duJi 



lui Jb—ioo 

Integrate e~^^^ over momentum space which is easy to do; this gives 



(2.23) 



n{E) 



(27rm) 



N 



Nl 



I rb+ioo ^I3{E-Hp) 



27ri Jb 



dp 



Ud^ 



(2.24) 



^Here we follow Padmanabhan 
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Then use Table 15.2 in Arfkin to calculate the Bromwich integral in the square brackets 
and replace Hp by its definition in ( |2.21| ): 



n{E)= ' -/ \E+-Y-, A Wd^Ti, (2.25) 

^ ' N\{lN-l)\jv\ 2^|r,-r,|j 1/ ^ > 

and finally change to new variables r = ri — r2 and with u,v = 2,3, ...,N in terms of which 

^ i2nm)-2^ r ^TT^3 (2.26) 
^ ^ iV!(|A^-l)!iy V" 

in this equality, 




Gm2 Gm2 1^ 



E + + V + - y r d^r. (2.27) 



3 iV^ 3 iV 

We now see that for r ^ 0, I — >oc / r^~~dr oc r^~~2". So for > 3 and r — > 0, $7 — > oo. 
Considerations of simple models show that the system with a subset of particles closely bound 
together by gravity and the remainder banging around with the high energy released has a 



large phase space volume associated with it [48]. Such a state can presumably be reached 



given sufficient time for the system to evolve. A system of point particles is however not 
realistic. In more realistic models, one make use of small hard spheres or fermions or remove 
the divergence of the potential with artificial cutoffs. This makes n{E) converge. 

It is important to realize that once converges by say taking hard spheres, the stable 
state is not necessarily a dense core with a dilute halo; other configurations may be more 
stable. For a sufficiently small cut off at least two particles will be very close together since 
such states would dominate were the hard spheres reduced to points. Other configurations 
may only be more stable if the cut offs are big enough. 



Models like these have been reanalyzed recently by Chavanis |12] who reviews earlier 
interesting works with various cutoffs. Neither of the models give however a direct evaluation 
of ^{E). 



(in) The meanfield approximation. 



A meanfield theory is a most common approximation when short distance effects are 
negligible on the timescales considered. A steepest descent calculation gives a meaningful 
statistical thermodynamics, equilibrium configurations and stability conditions. 

The following calculation is kept simple by not introducing any sort of cut off while 
embedding divergences into terms that would normally converge if there was a cutoff. These 
quantities when finite have no role in the results that interest us. 

One should however be aware that the maximum of the new integrand e'^ is not unique 
nor is it necessarily the highest. A second maximum may be much higher if the hard spheres 
are small enough. However, the calculation of that maximum is precisely brushed under the 
rug. A correct calculation is far too difficult. Our entropy will be associated with what is 
often referred to as a "local" maximum. This local maximum may not be the highest but it 
may be the more important one on time scales relevant in astronomy. 
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The calculation of i}(E) in the meanfield approximation is slightly intricate. We give 
here a heuristic derivation that will make the answer plausible. A complete and rigorous 



derivation that follows paper |25] is given in Appendix A. 

Let us ask how ( |2.10| ) should be modified when U{r) is replaced by the gravitational field 
of the particles themselves? The field is now variable with oo^ degrees of liberty. The correct 
O will have a functional integration over the space of all fields instead of a discrete integration 
over the rj-space. a needs a correction also that must depend on the field because the mean 



value of the energy in equation ( 2.11 ) is not the correct expression for a a self-gravitating 
system: the potential energy is counted twice. On the other hand the correction of a should 
not depend on a or /? since S must also be the correct expression for the Boltzmann entropy 
whether the field is a self-field or an external one; the local distribution does not know what 
the origin of the field is. 

With these remarks in mind, we now write what is the expression for the Gibbs density 
of states function in the "meanfield approximation" derived from ( 2.23| ) in Appendix A. 



First we must use a field variable W slightly different from and more convenient than U. We 
replace U by /3~^/^VF(r). In terms of VF(r) 

1 r+oo pa+ioo rb+ioo 

i^z,nij u J —oo J a—ioo Jb—ioo 

in this expression ;B is a constant which actually diverges if there is no short distant cutoff 
(see appendix A). VW is the element of volume in function space (see also appendix A) while 

ga-,3[^p2+m/3-V2i4.(r)]^^^ / ^FAH^dV . (2.29) 

SvrG J-oo 

A is the Laplacian. a differs from a in ( |2.1[1| ) mainly by the term WA.W. This addition pro- 
vides the correct mean value of the total energy as well as Poisson's equations for equilibrium 
configurations as we shall now see. And it does not depend on a, (5. 

(iv) A steepest descent evaluation. 

We start as in the case with an external field by calculating the extremal values of 
cr[a, /3, W^(r)]; these are associated with equilibrium configurations. 

— 1/2 

The equilibrium value We{r) of W{r) or rather U{r) = (5e ^e(r) is defined by 



e<^e-P.[^P +ml/(r)]^3^^ (2.30) 

-oo 



5W{t) 

p is the mean mass density of the particles in a meanfield U with a Boltzmann distribution 

/g = ffi^e-Pe[^P^+'mU(r)] ^ ^a^-f3,Eo^ (^2.31) 

Eq is here the energy of one particle in the meanfield U itself, a solution of Poisson's equation 
( 2.30| ). The equilibrium values {ue, Pe) of {a, (3) are fixed by equations similar to (2.11) but 



E is not quite the same as E: 

dcr ,^ /■ „ , , da 

— = 0^ N 



Jf.d. but i,=O^E = J 



da J dp 



1 1 
— p2 + -mU{r) 
2m 2 



fedco. (2.32) 
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The maximum cje of a, which is approximately the maximum of the integral when >> 1 is 
equal to the Boltzmann entropy. It is indeed easy to see that (|23l])(^) provide the same 
expression for the extremal value of a as ( ^.29 ): 



S(E, N,V) ae = - J fe Hfe)du + N. (2.33) 

We may thus write ( |2.28| ) as follows introducing 6^a which represents the sum of all the 
terms of order higher than one in a Taylor expansion near the extremum: 

1 r+oo ra+ioo rb+ioo 



oo J a— too Jb—ico 



Since, as before, the exponent is of order N the exponent in the integral is very steep and a 
good approximation will be obtained by limiting the integration to terms of order two. The 
reader interested in the derivation of d'^a and its explicit expression will find them both in 
Appendix B. There we show that ^{E) is approximately given by 

e^-'^^^^^-'^^^^^feBZ/g.) ^s'^'v{5W), (2.35) 

J —oo 

is a non-local quadratic functional in 5W which must be negative for the integral to 
converge, something of the form: 

/+00 
6WO{t, r')5W'd^rd^r' < 0. (2.36) 
-oo 

SW can in principle be expanded in terms of a complete set of orthonormal eigenfunctions 
say ^a(r) of the linear operator 0{r,r') with arbitrary discrete variables SY"". The quadratic 
form is then replaced by an infinite sum of squares {5Y"')'^ with coefficients that for 
convenience we write — ^A^; thus 



a=oo 

= -- E ^'aiSY^'- (2.37) 

^ a=l 

The A(j's can be arranged in increasing order X'l < X2 < X'^ < • • •■ Clearly the integral in ( p.35D 
will be convergent if X'l > 0; 5* is then a maximum and X'l > is a condition of stability. 
If >> 1 and 6"^ a' is sufficiently small, Sg c:t S as is usually expected from a steepest 
descent calculation. We can in principle evaluate the limit of validity of this quasi-equality 
by calculating 6^a' . 

(v) Convergence of Vl{E) in meanfield theory. 

There are two factors that cause il.{E) to be zero if it converges: B as we said diverges 
and the spectrum of eigenvalues is unbounded if the domain of definition of eigenfunctions 
is bounded. This is true under very general conditions that do not depend on the form 



of a finite domain It is true in particular for isothermal spheres [^] as shown in 

appendix B. In other volumes we may restrict 5W's to functions that satisfy a Poisson 
equation A5p = 4TrG5W and use these Jp's as variables; they are zero outside V. In that 
case, either > 0, the Boltzmann entropy is not a maximum of a which tends to plus 
infinity and S7(£') is undefined in this approximation because B in the denominator of ( ^.34)) 
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divergence. Or b'^a' < and the Gibbs density of states function tends to zero! There is no 
way out of this difficulty in a meanfield steepest descent calculation^. 

The culprit of the trouble is easily traced to the Laplacian operator whose inverse is the 
divergent 1/r potential Some form of short distance cutoff should have the effect of 

replacing the Laplacian by a non local operator with a convergent set of negative eigenvalues 

and a convergent B. 

Notice one important point: that the number of nodes of the eigenfunctions grows steadily 
as increases. When the spacing between the nodes in the finite domain is of the order of 
the short distance cutoff the Laplacian in 5^ a' is no more a reasonable representation for a 
convergent ^}{E). This means however that a short distance cutoff will changes the higher 
end of the spectrum and that X'l > is a perfectly good stability condition for the meanfield 
configurations with a Boltzmann distribution. 



(vi) Local and global maxima. 



Had we made a correct calculation with some specific cutoff not too big, we would most 
likely have found at least two maxima of a, the one we have just found and another one 
associated with a highly concentrated core and entropy 5*0 say. This Sq would tend to 
infinity if the cutoff went to zero. A steepest descent calculation would have given 

n{E) = e^o ~ e^«+^ (2.38) 

Since both S and Sq are of order >> 1 the Gibbs entropy, for most values of the energy 
E and volume V would be approximately equal to either S or Sq. Such a transition appears 
in a Monte-Carlo calculation of i^{E) by de Vega and Sanchez [17|. 



In principle the global maximum is the equilibrium configuration after an infinite time. 
However, the system may find itself in the local maximum for a very long time indeed. This 
makes local maxima interesting. 



(vii) Remarks about relativistic systems and other ensembles. 



(a) Relativistic systems. 



The Gibbs density of states function can and has been calculated in general relativity 
The advantage of general relativity is that it is a mean field theory to start with. 
The disadvantage is that the formalism is far more complicated, arbitrariness of coordinates 
creates additional difficulties in functional integrations and moreover a density of energy is 
not defined. With so much trouble, one finds however that the total entropy is again given 



by ( 2.33 ) but of course fe is different from ( 2.3f 



Other ensembles have not found wide applications in astronomy but they are of theoretical 
interest. 



(b) The grand canonical ensemble. 



In this ensemble N and E are not kept fixed. What is fixed is a and /3. It follows 
that the Gibbs equivalent for this ensemble I^gCi called the grand partition function, is in 

®see also the end of Appendix B in this respect. 



9 



the meanfield approximation given by expression ( ^.28| ) in which da and ^e^^dp 

integrations are removed from e'^ because and E are not fixed anymore; thus 

1 f+oo 

^GC = / e^'G^VW =^ (JGC = + aN - (5E. (2.39) 

flac is thus a double Laplace transform of Q{E). Thermodynamic equilibrium is associated 
with the extremum acce of the exponent ctgc- This extremum is simply related to the 
entropy: 

acce = S + aN - I3E. (2.40) 

The thermodynamic potential has no conventional name. N and E are now mean values 
which are however related to a, /3 by the same equations ( 2.32| ). In this ensemble like in 
the following ones equations ( p.30|) , ( p.31| ) and ( p. 32 ) always hold at the point where the 
exponent is extremum. But equations ( |2.32 ) have slightly different interpretations. 

(c) The canonical ensemble. 

In this ensemble the system has a fixed number of particles and the temperature or (5 is 
kept fixed rather than the total energy. The Gibbs equivalent for this ensemble is called 
the partition function and it is obtained from ( 2.28 ) by removing the ^G.^^df3 integration 
from e'^ . Thus 

1 ra+ioo r+oo 

Qc = — ^ / / e'^'^daVW ac = (J-(3E. (2.41) 



2mB 



oo 



The corresponding thermodynamic equilibrium function which plays the role of 5 and asso- 
ciated with the extremum ace of the exponent ac is 

ace = S-f3E = -(5F. (2.42) 

The thermodynamic potential F is known as the free energy. 

(d) The grand micro canonical ensemble. 

This fourth ensemble not often considered |44| is one in which the total energy is fixed 
but the particle number is allowed to fluctuate. This ensemble is physically plausible if there 
are pair creations of particles. Thus a and the total energy are fixed. The Gibbs equivalent 
^GMC is obtained by removing the da integration from e°": 

1 rb+ioo r+oo 

^GMC = ir-E / e'^^^'^dpVW aGMC = <y + aN. (2.43) 

I-KIO Jb-ioo J-oo 

Other ensembles have been studied like isothermal spheres under constant pressure. 

(viii) Summary and conclusions. 

Here now are the main points of this section. 
(1) We started from the Gibbs density of states function Q, equation ( |2.22| ), which we wrote 
in a meanfield approximation as equations ( 2.28| ) ( p. 291 ). 
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(2) We used a steepest descent method to calculate equilibrium configurations which are given 
by Poisson's equation ( |2.30| ) with a Boltzmann distribution of particles in phase space (2.31). 
We also calculated the entropy ( 2.33| ) and got a meaningful expression if ( 2.34 ) is convergent. 
The condition of convergence is a condition of stability. The results are only valid if some 
sort of short distance cutoff is assumed, the details of which are however unimportant. The 
limit of validity of the whole approximation scheme can be evaluated if V is specified. It is 
not clear however if the Boltzmann entropy is a global or a local maximum. That depends 
on the cutoff and it would show up in a correct calculation. 

Lynden-Bell and Wood, in dealing with isothermal spheres proceeded in a simpler way. 
They started from the entropy defined by ( |2.33| ) with an unknown distribution function / 
restricted by conditions (2.32) in which j3~^/'^W = C/ is a solution of Poisson's equation 
(|2.30| ) with p = m J fdu). This out of equilibrium entropy is then extremized {5S = 0) with 
arbitrary (5/'s subject to two restrictions 5N = 5E = 0. The extremum of / is the Boltzmann 
distribution function ( 2.31 ). Padmanabhan (5^] used from the same method to find stability 
conditions. This method is simpler and well adapted to gaseous systems. It short-circuits 
many of the difficulties encountered above. Thus the statistical thermodynamics of Gibbs 
ends up with the same mathematics but it must be said again that "it has a particular beauty 
of its own [and] is applicable quite generally to every physical system" [30|. 

Section 1 dealt with the thermodynamic equilibrium of N point particles of mass m 
isolated in a volume V with energy E and entropy S. Section 2 deals with thermodynamic 
functions which give statistical thermodynamics its physical content. We also develop the 
thermodynamic stability theory of equilibrium configurations. The theory of stability in 
section 2 applies to a much wider class of systems than the one studied in this section. 



3 Thermodynamic equilibrium and stability 

The first subsection deals with the thermodynamics of the system of N particles with grav- 
itational interactions in a volume V . This is our generic case. The rest of section 2 is more 
general. Please notice that we shall be concerned by equilibrium configurations only; we drop 
therefore the index e from Og, (3e and /e which is we have no reason to use anymore and write 
simply a, (5 and /. 

(i) Thermodynamic functions of self gravitating particles. 

Besides E, N, V and S{E, N, V) there are other thermodynamic functions that give a 
physical content to the results. Consider in particular the derivatives of S with respect to 
E,N,V. The derivatives are easily calculated by considering a defined in ( |2.29| ). This is a 
function of E, N, V as well as of a, (3, W{y). But the entropy is an extremum of a for which 
da I da = dajdfi = 6a/6W{r) = 0. Therefore the partial derivatives of S with respect to 
E, N, V are the same as the partial derivatives of a with respect to E, N, V keeping a, /?, W{r) 
fixed. The partial derivatives must of course be evaluated at the extremum. A quick look at 
( 2.2g| ) shows immediately that: 



dS ^ dS ^ dS f f+°° 3 \ pb 
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the index means "on the boundary"[]. It is worth noticing that the local pressure 

^ _™V^V = ^ /«V = thus - = m. (3.2) 

where Pi, is the pressure on the boundary. It follows from ( |3.lD and ( |3.2| ) that 

1 a 

dS = f3dE - adN + pPbdV or -dS = dE - -dN + PbdV. (3.3) 

This differential expression, the most important one in statistical thermodynamics, shows 
that 

1 a 

— = T is the temperature, and — = /i is the Gibbs chemical potential. (3.4) 

P P 

T,n.Ph are global thermodynamic quantities; they concern the whole system. 

(ii) Conjugate thermodynamic functions. 

Of particular interest in stability analysis are pairs of conjugate thermodynamic functions 
with respect to the entropy like {E,f3 = dS/dE), {N, -a = dS/dN) and {V,(3Pb = dS/dV). 
Pairs of conjugate parameters in other ensembles are easily obtained from ( |3.3| ). For instance 
in the grand canonical ensemble the pairs of conjugate parameters with respect to ctqcc are 



{I3,~E), {a,N) and {V,(3Pb) since from (|J) and (MQ) we see that 

daoce = -Ed(3 + Nda + pPbdV. (3.5) 

In the canonical ensemble pairs of conjugate parameters associated with the free energy or 
rather —^F are (/?, —E) and (iV, —a) and {V, f3Pb) and for the grand micro canonical ensemble 
the pairs of conjugate thermodynamic functions with respect to acMCe are {E,j3), {a,N), 
{V,PPb). From here on the thermodynamics of self-gravitating systems can be developed 



along the royal path taken by Lynden-Bell and Wood |48] for studying isothermal spheres. 
Now we continue at a more general level. 

(Hi) Stability conditions for equilibrium states in general. 

If the Gibbs density of states function in the mean field approximation is convergent and 
the steepest descent evaluation of Q{E) to order two is meaningful, equation ( p. 34 ) together 



with equation ( 2.37 ) suggests that with a proper choice of variables in 0(£'), say x"", in ( ^.28 ) 



the stability conditions would appear directly in their simplest form. Thus starting from 
any Gibbs density of states function of any ensemble in these appropriate variables we would 
write ^}{E) in a form like this: 

n{E) = i /" e'"^^'^") (3.6) 

where C is a constant that depends on the choice of the x'^'s. The extremum of w which 
defines equilibrium configurations is then of the form 

dw 

— = whose solutions = X^iE); (3.7) 



^ There would be no mechanical equilibrium otherwise as the boundary must be an equipotential. Indeed, 
Pb must be t he s ame everywhere on the boundary as mechanical equilibrium requires: —VP — pVU = which 



implies, see (3^), U (xlnp and Ue must be the same everywhere on the boundary because the pressure force 



is normal to the surface. 
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there may be more than one solution for given E and whatever other parameters there may 
be, say C,q, {q = 1, 2, ...). Near an extremum where x"" = X"'{E), the expansion of w would be 
of the following form: 

w = We-^J2^-i^"-^''i^)? + 0s. (3.8) 



Equations (3.7) and (|3.8D are similar to those obtained for testing equilibrium and stability of 



a mechanical system with potential energy —w, variables and parameters {E,C,q) though 
the number of variables in this case is usually finite. The whole analysis as can be seen is 
in fact valid for self-gravitating systems or other systems with short range forces in thermal 
equilibrium. 

The calculation of the stability limit Ai is not easy in general. Therefore the following 
criteria may be helpful. 

(iv) A general criterion of stability. 

(a) A useful identity 

For definiteness we keep in mind the micro canonical ensemble and the w{= a) intro- 
duced in the previous section but keep our new notations that refer to more general systems. 
Consider the parameter^ E and introduce its conjugate with respect to w rather than with 
respect to We- This defines a sort of inverse "temperature out of equilibrium" say /3(£';x"). 
The derivative of w given by equation ( p.^ ) with respect to E keeping all x'^'s and Q's fixed: 

P{E;x^) = — = P{E)+J2Xa^[x'^-X^{E)]+02 , m = q^- (3-9) 

Clearly the equilibrium value of /3 is /3 and the second order terms is to be interpreted as 
fluctuations of /? near equilibrium. When w = a, (3 is the inverse temperature; here (3 is some 
"generalized" inverse temperature. 

Let us calculate first order derivatives of f3 with respect to x"" and E at the point of 
extremum. At that point derivatives of O2 in (^]^) are at least of order 1 and will thus add 
nothing to the derivatives at the point of extremum. The second terms of (3.9) contribute 



only through derivatives of [x"" — X'^{E)] because this quantity is zero at the extremum. One 
thus obtain: 

= Aa { „ \ no summation on a , ( 7-^ 1 = 7-^ — \ J Aq ( ) . (3.10) 





V dE 

Extracting dX^/dE from the first equality and inserting into the second gives the following 
expression for the slope of the "linear series" of conjugate parameters P{E) with respect to 

We- 

df3 _ (df3\ {dl/d^ 

dE-\dE)^^\^^^- ^'-''^ 

Let us emphasize again that the calculation is valid for any pair of conjugate parame- 
ters with respect to any w, being the statistical weight of any thermodynamic system, 
including one with short range forces between the particles. It also applies to any pair of 



'Any other parameter might do. 
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conjugate parameters with respect to (minus) the potential energy of a mechanical system in 
equilibrium. It is however extremely useful to keep in mind the concrete example of the 
gravitating particles of mass m with energy E in a volume V. 



(b) Poincare's criteria of stability. 

There exists a classical result of Poincare [^ ] which is discussed in various treatises on 
stability, see for instance Jeans p^] Ledoux ||45(| or Lyttleton E^, and was applied to self- 



gravitating systems in thermal equilibrium for the first time by Lynden-Bell and Wood |48]. 
It says the following. A change of stability may only occur where two or more linear series of 
equilibria like P{E) for instance]^ have one equilibrium configuration in common ("bifurcation 
points") or where two or more series merge into each other ("turning points"). When this 
happens stable equilibria may turn into unstable ones; reciprocally, unstable ones may become 
either stable or more unstable. Thus a change of sign of Ai and perhaps of more eigenvalues 
or stability coefficients appear only at bifurcations and turning points. 

The Poincare method of linear series does not say whether a change of stability actually 
occurs or not. It has been pointed out |]3l|| that under very general conditions the stability 
or the number of unstable modes (the number of negative Aa's) can be deduced from the 
topological properties of series of equilibria, i.e. from purely thermodynamic considerations 
as we shall now review. 

(c) Simple eigenvalues. 

Suppose the spectrum of eigenvalues is simple that is non-degenerate: Ai < A2 < A3 

This assumption is far more than we shall ever need but it makes further explanations simpler. 
In practice only some of the smaller Aa's need to be distinct. Consider a line /3{E) with 
points that correspond to stable configurations (Ai > 0) and follow the line towards a limit 



of stability. As we approach that limit Ai — > the 1/Ai term in (3.11) begins to dominate 
and 

II ^ {dp/dx ), ^ provided {dp/dx^)e^O. (3.12) 

oE Ai 

We shall come back later to what happens if {d$/dx^)e —>■ 0. If equation ( |3.12| ) holds, the 
slope of (3{E) has the same sign as Ai. As a result: 

(1) Two linear series of conjugate thermodynamic functions merge into each other at a 
point where there is a vertical tangent. A change of stability occurs thus necessarily there, 
at a turning point. 

(2) The stable branch is the one with a positive slope near the vertical tangent. Another 
way to state this property is as follows. Stability is lost where the branch turns counterclock- 
wise. A stable branch never turns clockwise. 

(3) Say Ai < and A2 ^ along the P{E) line. There is thus already one mode of 
instability. If a further mode of instability shows up it must be at a new turning point where 



A2 — > 0. Near the new turning point the dominant contribution to ( 3.11 ) comes from the 
I/A2 term and equation ( p. 121 ) holds with an indice 2 instead of an indice 1. 

(4) Say Ai < < A2 and P{E) turns clockwise. This means that Ai ^ again and that 
unstable configurations turn into stable ones on the side with positive slopes. 



'Any linear series, not only of conjugate thermodynamic functions like P{E). For instance, one of the X°''s 
as a function of one of the (q^s or a P{V) diagram. 
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(5) If a stable linear series spirals inwards counterclockwise, stability is lost at the first 
turning point. We give an example of this case in section 4. 

(6) If a change of stability occurs at a point where {d/3/dx^)e — > and Ai ^ the 
previous considerations become obviously invalid. Notice however that a small perturbation 
of the potential energy easily removes the coincidence. Indeed, let us change w to w + eEx^ 
say with |e| << 1. The equilibrium equations ( |3.7D become 




dw _ , , , dw , , , , 

— + eE = {\e\«l) and _ = (a > 1). (3.13) 

while 

dEdx^ ] ^ dE Ai ^ ^ 

/ e 

This shows that such coincidences will show up in the mathematics but not in experimental 
observations [^] nor are they likely to appear in numerical calculations and many interesting 
models in astronomy and astrophysics are solved numerically. 

(7) Finally if we know that one configuration is stable or we know its degree of instability 
we shall know about the stability or degree of instability of every other configuration of the 
linear series. 



(d) Bifurcations and degenerate spectra. 



Bifurcations like degeneracies of the spectrum of stability coefficients come often from 
an excess of symmetry. A characteristic property of bifurcations in problems of physical 
interest is that they are "unstable". Small perturbations of the potential energy as just 
indicated remove unstable bifurcations and transform them into turning points [^3|; such 
bifurcations do not show in experiments or in numerical solutions. We shall give an example 
of an unstable bifurcation in section 4. The theory of bifurcations^ is of some interest in 
mathematical physics Q. 

In the rest of this section we deal only with simple spectra and turning points at vertical 
tangents. 



(v) Cy < and the nonequivalence of ensembles. 



(a) About Cy < 0. 



We can now see that in a thermally stable isolated system but near instability (Ai 0) 
the heat capacity Cy = dE/dT is negative in stable configurations and turns positive in 
unstable configurations by going through zero. Indeed, according to ( 3.12| ), 



^"The problem of distinguishing stable from unstable bifurcation offers little difficulty with linear series 
depending on one parameter j^. The unfloding of bifurcations with more than one parameter and several 
variables is much more complicated. Thom classified all structurally stable inequivalent unfoldings or as 
it is called "elementary catastrophes" with up to two variables and four parameters. Great progress was done 
in this subject since 1975. 

Catastrophe theory has not been very useful in physics. It did not contribute a thing that was not already 
known otherwise. It provided however a sound mathematical basis to the theory of bifurcations and to the 
classification of inequivalent ones. 
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The same will of course happen in an unstable system on the verge of acquiring a second 
mode of instability (Ai < < A2 ^ 0). Cy may thus also be negative in unstable systems. 
The sign of Cy is therefore not a criterion of stability like in classical thermodynamics. 

The fundamental reason for this difference is the fact that the system is no subsystem 
of an ensemble. This is made obvious by the following standard explanation of why small 
subsystems are stable if their heat capacity is positive. 

Consider an isolated system of energy E made up of two systems of energy Ei and 
E2 = E — El of comparable size. Let Si{Ei) and 6*2 (£'2) be their entropies. The sum of the 
entropies is the entropy of the system out of equilibrium but for simplicity we write 



SiE)=SiiEi) + S2iE2). 



(3.16) 



Equilibrium is reached when 
-fdSi 



6S 



\dEi 



fdS2\ 

\dE2J, 



5Ei = (/3i - f32)6Ei 







Pl=P2 = p. 



(3.17) 



To check stability we calculate second order variations of S which at equilibrium must be 
negative: 



dEi 



+ 



dfh 
dEn 



{SEi 



1 



(3^ 

2 \Cvi 



1 1 

+ 



Cv2 



{5Eif < 0. (3.18) 



Since SEi = CyiSTi = —CyiSPi/P'^, we can rewrite the last expression as follows: 



dfh 



< 0. 



Now if system 1 is a small subsystem, Ei <^ E2, also Cyi <C Cv2 and 

2 



(^'S)e^-\Cyi(^-f) <0. 



(3.19) 



(3.20) 



The stability conditions ( p. 20 ) implies that a stable small sub-system must have a positive 
heat capacity as is well known. 



(b) About the non-equivalence of ensembles. 



If instead of keeping the system isolated we keep it at constant temperature, the role of S 
is then played by the thermodynamic potential wce or —f3F, see equation ( p.42| ), and (/3, —E) 
is a pair of conjugate parameters with respect to —f5F. But —E{f3) is the same line as f3{E) 
drawn in coordinates rotated 90° clockwise. Vertical tangents become horizontal ones and 
reciprocally. Let \ac be the Poincare coefficients of stability for the canonical ensemble. Near 
instability Aic — > 0. The analog of ( 3.12| ) is now 



d{-E) 

df3 



d{-E)/dxl; 



Aic 



This shows that near instability a stable canonical ensemble has Cy > 0. Nevertheless 
unstable configurations may also have Cy > 0. For instance if Aic < < A2C and A2C — * 0. 
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Thus Cy > is not a criterion of stability for a canonical ensemble either. It is surely a 
sufficient condition of stability near and indeed far away from instability. 

The canonical ensemble and the isolated system or microcanonical ensemble cannot be- 
come unstable for the same equilibrium configuration unless vertical and horizontal tangents 
appear at the same point. This difference between ensembles is referred to as the non- 
equivalence of ensembles and is a property of non-extensive systems like self-gravitating ones. 
In extensive systems, different ensembles are equivalent Why does the canonical ensem- 
ble get unstable and the other does not is nicely explained in the Lynden-Bell and Wood 
paper [p^]. 



(vi) Summary and conclusions. 



Under general conditions, mainly that the spectrum of Poincare coefficients of stability 
be simple, stability limits and the number of unstable modes can be found for all calculated 
configurations using only linear series of thermodynamic functions of equilibrium configura- 
tions provided we know whether one configuration is stable or what is its degree of instability. 
Conjugate thermodynamic functions are particularly valuable in this respect as changes of 
stability show up at vertical tangents and the slopes in their vicinity have simple interpreta- 
tions. 

The thermodynamic criterion has a number of limitations. One rarely calculate all the 
sequences of equilibrium and therefore some bifurcations may not show up because the branch 
points are missing. Thus equilibria might become unstable and the system might choose to 
be in a more stable state which has not been calculated. A complete answer needs a detailed 
analysis of the second order terms of the Gibbs density of states function like the one started 
in subsection 2{i). The much simpler version of Padmanabhan |54] is also useful. 

A less severe limitation is that the method is mathematically not a full proof because we 
must assume that the spectrum of stability coefficients is simple and at the same time that 
conditions like ( ^.12 ) hold. The method may however be applied to numerical solutions with 
a reasonable degree of confidence because this type of mathematical singularity will not show 
up in numerical calculations nor indeed in experiments. 

Another general limitation of any thermodynamic criterion of stability is that we learn 
little about the nature of instabilities, triggering mechanisms, and what becomes of stable 
states which evolve through a series of quasi-equilibria along the linear series up to and 
beyond the limit of instability. 

We have seen that the sign of heat capacities is not a criterion of stability and that in 
microcanonical and canonical ensembles near instability heat capacities have opposite signs. 
These ensembles are not equivalent as in classical thermodynamics. The proofs of these 
properties have a great degree of generality. The next section deals with fluctuations. 



4 Fluctuations 

(i) Fluctuations in self- gravitating systems. 

The thermodynamics of fluctuations in stable equilibrium configurations becomes inter- 
esting near instability. Fluctuations can put the system out of equilibrium. Calculations of 
fluctuations are particularly interesting when dealing with metastable states of systems like 
the system of point particles in an enclosure considered in section 1. The entropy of such a 
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system may be a "local" maximum smaller than the absolute maximum. In that case there 
are two maxima and between them there must a saddle point. Such a saddle point must show 
up close to a limit of stability because there the linear series turns counterclockwise. Figure 
1 illustrates the situation with a linear series P{E) for isolated systems. Let /3 be the inverse 
temperature of a stable system and j3' < (3 the inverse temperature of the unstable system 
with the same energy E. If the mean quadratic fluctuation of temperature or rather of its 
inverse < > is equal to or greater than the square of the gap (/3 — the statistical 

weight of a real fluctuation bigger than (/3 — [3') becomes important; this may put the system 
in a state where w can increase to a bigger local maximum or to the global maximum. That 
is the entropy will have a chance to increase which it will surely do and equilibrium will be 
lost. It is thus interesting to evaluate the statistical weight of the fluctuations. Fluctuations 
can only displace the limit of stability towards higher values of (3 and the question is how 
much higher. This is the subject of this section. The theory which has the same degree of 
generality as the theory of stability of section 2, comes from works by Okamoto, Parentani 
and myself [53| p6|. 



(ii) The probability of a fluctuation near instability. 

Figure 1 is a good starting point. The point of marginal stability has coordinates {Em, Pm) 
and a vertical tangent. Consider near that point a stable configuration with coordinates 
{E,/3). The dominant contribution to Q near instability, see (p.6|), must come from 6x^ in w 
since — > cxd when Ai — > 0. Let us therefore integrate over all variables except x^: 

-I r+oo a=oo s a=oo \ 1 p+oo 

n(E) ^s^kT.MS-"? Yldx-^^-Xir-f) e~-2^^^'^'y'di5x'). (4.1) 

As long as the exponent is sufficiently steep (thus not too close to where Ai = 0) we may 
integrate over 6x^ between ±00 as indicated. 

Notice that Sx^ = x^ — where is the equilibrium value of x^ at point {E, /3) and x^ 
itself is some value out of equilibrium with the same energy E. We have a relation between 
6x^ and the fluctuation of temperature 6P = /3(x^; E) — (3 which to first order is simply 



5/3 = /3-/3~ 5x^. (4.2) 




Replacing {dj3/dx^)e in terms of dE/dp using the second of equations (|3.12D , 

,,^^^JL^^U^?EV'\^.n^ (4.3) 



We use both these equalities to replace 6x^ in the integral of (4.1), the first equality in d{6x^ 



the second in the exponent, and obtain 



We conclude that the statistical weight of a fluctuation of f3 to between /? and (3 + dp is 
proportional to the integrand of ( [4. 41) . Thus the probability dV of the fluctuation is given by 
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the normalized expression: 



dV = ^e-''dt with t' = l^(p-pf. (4.5) 
^/7^ 2 op 



The mean quadratic fluctuation of temperature is thus given by 

Here we used ( p. 21 ) to get the expression on the right hand side. The formula is like that of 
Landau and Lifshitz |Q for mean quadratic fluctuations of temperature in a small subsystem; 
it has by necessity the opposite sign for Cy- 

Fluctuations of other mean values near instability in different ensembles can be derived 
in a similar way. 

The formula ( |4.5| ) may be made more specific because we are near the turning point 
where {dE/df3)m = 0. {dE/d(5) is thus small and to the lowest order in (/? — /3m), assuming 

{d^E/dp^)m + 0, 



t^ defined in ( [4.5D can also be written 



/3"^ 



(4.7) 



,2 ^(rr.QCy\ ff3-f3m\(P-P 



2 



and the mean quadratic fluctuation satisfy the following equation: 

If {d'^E / dpi^)m = a higher even derivative must be different from zero because £^ is a 
minimum and the final results will not be very different; there will be higher powers of 

fi-Pm. 

(Hi) Stability limits induced by fluctuations. 

We now define a new limit of stability, the point {(3i,Ei) on the linear series in figure 1 
where the mean quadratic fluctuation equals the square of the difference between (3 and (3' 
for reasons explained at the beginning of this section. Near the point of marginal stability 
{(3m, Em) we have approximately {f3i — f3[) ~ 2{f3i — 13m)'-, thus f3i is defined by the condition 
that 

<0-f3i)^>^A{(3i-f3mf. (4.10) 



We shall obtain f3i by inserting this value of the mean square fluctuations into (|4.9|) ; this 
gives 



dCv\ fPi-P^, 



3 



(4.11) 
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Notice that Cy is of order N and therefore the change from fim to (3i is of order N~^/'^ which 
is not necessarily very small. The corresponding value of Ei is readily found by expanding 
E{j3) in a Taylor series to order two near f5m'- 

To have some idea on how sharply defined the point {(5i,Ei) is, consider an equilibrium 
configuration with [3 > [ii and evaluate the weight of the probability for a fluctuation to 
induce instability, i.e. calculate e~* for 

/3-/3~2(/3-/3^). (4.13) 

If we insert this value for (3 — (3 into ( |4.8| ) , taking account of ( [4.11 ) we find that 

,2 l^/5-/3„ " 



For P — Pm = Pi — Pm, e~*^ ^ 0.6 but if P — Pm is twice as big as Pi — Pm-, e"*^ ^ 0.02 and 
if equal to three times that difference, e ~ 10~ . The exponential fall is relatively sharp 
and the new limit of stability (Ei,Pi) induced by fluctuations is not too badly localized. 



(iv) Summary and comments. 



Fluctuations in non-extensive systems like self-gravitating ones have a peculiar behavior 
which has no analogue in classical thermodynamics where smooth turning points do not exist 
and ensembles are equivalent. Fluctuations can induce instability at lower temperatures 
(higher /3's) than the theoretical limit Pm- The relative change is of order N^'^^'^. N = 10^ 
is often considered in numerical calculations to study globular cluster models. Instability 
limits found in these models may be quite different from more realistic ones with A'^ = 10^ 
or higher. Moreover, since A^~^/^ is not necessarily a big number the stability limit may be 
sensitive to the values of {ATdCy / dT)m- 

Fluctuation theory is useful when slow evolution occurs towards instability as we shall 
see in section 5. 



5 Examples 

In this section we give two examples in which thermodynamics of self-gravitating systems is 
used. The first one is that of isothermal spheres. Equilibrium configurations were obtained 
numerically. In this case linear series of conjugate thermodynamic functions correctly give 
stability limits. The calculation of fluctuations give useful additional information. The 
second example is that of liquid ellipsoids. Here equilibrium configurations are defined by 
quadratures and the multiple linear series are replete with bifurcations. The bifurcations are 
however unstable and may be lifted with small perturbations of the potential energy. Other 
examples are mentioned in less detail. In some of these examples, the stability theory of 
section 3 has been used and in others it has not but might have saved complicated calculations 
to obtain the same results had it been used. 

(i) Isothermal spheres. 
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Isothermal spheres belong to standard literature in astronomy [Q]. This well posed math- 
ematical problem can be solved numerically and provides a not too crude model for cores of 
globular clusters and played an important role in their understanding. 

(a) Equilibrium configurations. 

The problem at hand is the one treated in section 2(iii) with a spherical enclosure of 
radius say. Antonov [Q] has given an elegant proof that in a sphere, spherical configurations 
maximize the entropy. Spherical equilibria have a gravitational potential U which depends on 
the distance to the center r only and Poisson's equation ( |2.30| ) with a Boltzmann distribution 



( 2.31 ) reduces in this case to 

/27rm\3/2 



dU{r) 



dr 

The boundary conditions for regular solutions are: 



dU\ , , GNm -GM , , 

— =0 and \J(rB) = \JB = = • (5.2) 

dr ) r-B TB 



Equation (|5.l| ) has a singular solution 

1 

27rGm/3r 



(5.3) 



Regular solutions oscillate around this singular solution to which they become closer and 
closer as r grows. 



Excellent tables of numerical solutions have been given by Emden |15], more detailed ones 
are found in Chandrasekhar |^ with some additions in Chandrasekhar and Wares [^. The 
tables are still useful today even if Mathematica provides 10^ exact digits in a fraction of a 
second for solutions with density contrasts up to 10^ in no time. 

The density contrast TZ = p{r = 0)/p{r = tb) is a useful parameter which characterizes 
the relative strength of gravity. It is obtained by direct integration of equation ( [5.1D because 
[/ oc Inp. Equilibrium configurations have been studied in detail in Chandrasekhar Q but 
also in Lynden-Bell and Wood [48| and in Padmanabhan |55]. A short review on isothermal 



spheres relevant to astronomy is given in Binney and Tremaine 

(b) Stable configurations. 

(q) Isolated spheres and isothermal spheres in a heat bath. 

Antonov [Q] found analytically that stable isolated systems can exist for density con- 
trastsQ 7^ < 709. A nice and simple analytical proof has been given by Padmanabhan |p 



Lynden-Bell and Wood [48| found the same result using a combination of Poincare's turning 
point method and direct calculation of the entropy. They also found that equilibrium con- 
figurations maintained at constant temperature are stable for IZ < 32.1. It is worth noting 
that Padmanabhan's method allows to calculate the density perturbations 5p that trigger 
instability. This gives interesting additional information that is beyond thermodynamic's 
capability. 

^^In general we give numerical results with 3 digits. 
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The method of hnear series with conjugate thermodynamic functions give these results in 
one stroke and adds new interesting information. Figure 2 represents /3m|[/B| as a function 
of E/M\Ub\] N, rs and thus Ub are constant. 

For small /3's, i.e. high temperatures, gravity plays a minor role, the particles behave 
like a perfect gas and perfect gases are stable for arbitrary perturbations [^]. The curve 
of figure 2 spirals inwards counterclockwise towards a point with coordinates E/M\Ub\ = 
— J, /3m|C/B| = 2. If we follow the stable branch from low density contrasts, say from E = 
where TZ ~ 6.85, to high density contrasts and negative energies, a vertical tangent appears 
at 7^ ~ 709. Thus isolated isothermal spheres become unstable for TZ > 709. This confirms 
once again Antonov's result but figure 2 also shows that all configurations with a density 
contrast TZ > 709 are unstable with respect to spherically symmetric perturbations and that 
the number of unstable modes increases with TZ because the linear series spirals inwards 
counterclockwise. 

Notice as expected that stable configurations near the vertical tangent at 7^ ~ 709 have 
negative heat capacities. Notice however that unstable configurations with 5.22 • 10^ < TZ < 
4.50 • 10^ have also Cy < though they are unstable and have one mode of instability: 
Ai < < A2. 

Notice finally that all configurations are stable for non-spherically symmetric perturba- 
tions 1 25 1 but this has not and can not be proven by thermodynamic arguments since we have 
linear series of spherically symmetric equilibria only. 

If we turn figure 2 clockwise 90° we are looking at —E{f3). This is the appropriate linear 
series for canonical ensembles as we have seen in section 2. We can follow again the line 
of stable configurations from E = Q and TZ = 6.85 towards higher density contrasts with 
negative energies, we shall meet a vertical tangent at TZ = 32.1. Thus isothermal spheres in a 
heat bath become unstable for TZ > 32.1. This confirms Lynden-Bell and Wood's result and 
shows in addition that all configurations with higher density contrasts are not only unstable 
but more and more so as —E{P) spirals inwards counterclockwise. 

Notice as expected that stable configurations near the vertical tangent have Cy > but 
unstable configurations with two modes of instability Ai < A2 < < A3 and density contrasts 
4.50 • 10^ < 7^ < 5.45 • 10^ have also Cy > 0. 

Chavanis used a method similar to that of Padmanabhan to calculate the density pertur- 
bations 6p that trigger instability in isothermal spheres at constant temperature [13| as well 
as in the ensembles studied in the next subsection [15|. 



(/3) Grand canonical ensembles and grand micro canonical ensembles. 



The method of linear series with conjugate thermodynamic functions has also been applied 
to systems in which neither energy nor the number of particles are fixed. In grand canonical 
ensembles a, /3, rs are fixed. An appropriate linear series for this ensemble is, see section 
3, the curve —E{j3) at fixed a and r^. This linear series, see also [^], is shown in figure 
3A and 3B and must be looked at rotated 90° clockwise. What is represented is actually 
(5* = P/(3o as a function of E* = {(3^Gm^/rB)E with /3o = 27r3m7(87rGr|e°)2. For S > , 
the non-dimensional quantities are of comparable magnitude. Notice that near E = there 
is a cross over but not a bifurcation with two different solutions whose stability or instability 
is known. For E < the curve winds in counterclockwise staying very near the vertical axis. 
The limit point has coordinates {E* = —5.37 • W~^, f3* = 13.6). This is why the diagram has 
a different scale for E < 0, figure 3a, and E > 0, figure 3b. 
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At low density contrasts and high temperature the energy tends to zero; kinetic and 
potential energies are both high and of comparable magnitude. It has been shown |24] by 
calculating fluctuations oiwcc near the extremum, see ( 2.40 ), that grand canonical ensembles 
are stable with respect to arbitrary perturbations^ for TZ < 1.58. We can see on figure 3 that 
grand canonical ensembles have indeed a turning point at 7^ = 1.58, —E{(3) has a vertical 
tangent there and turns counterclockwise. At the next vertical tangent IZ ~ 106 and the 
curve turns again counterclockwise and at higher density contrasts the curve is spiralling 
inwards. There is thus no stable equilibrium configuration for IZ > 1.58. 

The situation is more eventful in the grand microcanonical ensemble in which E, a, 
are fixed. Here figure 3a and 3b are appropriate linear series for checking stability. At low 
density contrasts equilibrium configurations are certainly stable: the ensembles are more 
constrained than the grand canonical one which is stable. We see that stable configurations 
exist for density contrasts IZ < 1.66. Configurations with 1.66 < TZ < 11.6 have one mode 
of instability. There exists a second series of equilibrium configurations with 11.6 < TZ < 
92.6. At higher density contrasts the curve spirals inwards counterclockwise and equilibrium 
configurations become more and more unstable. 



(c) Fluctuations. 



Fluctuations displace the limits of stability towards lower density contrasts. For instance 
in isolated isothermal spheres it was found that quadratic fluctuations would induce 
instability at density contrasts TZi ~ 709 • e"'^ '^'''^ . If N = 10^, a number close to 
numerical experiments, TZi ~ 510. The formula is still marginally valid for A'^ = 10 for which 
TZi ~ 154. These are not small effects. 

In this connection, it is worthwhile recalling Monaghan's [^] application of the theory 
of hydrodynamic fluctuations to a self- gravitating gas. He showed that density fluctuations 
become large before the point of ordinary stability is reached. Thermodynamics supports 
Monaghan's finding. 



(d) Relevance to astronomy. 



The simple model of isothermal spheres played an important role in understanding the 
structure of globular clusters. The model taken seriously by observers |5C] is that of Michie 
|51| which was put to extensive use by King |^^. The Michie-King model is a truncated 
Boltzmann distribution which in our notation is this, see ( 2.31J ): 



/ = A{e^^° - e^^=) for Eq < E^ and / = for Eq > E^ 



(5.4) 



Ec is an energy cutoff that simulates the absence of high energy escapers. Stability lim- 



its for this model and variants thereof |32| do not change the general trend of equilibrium 
configurations |35]. Linear series are always counterclockwise inwinding spirals and stable 



configurations exist up to some maximum density contrast that varies not very much. 

The distribution (^.41 ) is a simple model for classifying observable parameters of globular 
clusters. 

It is not clear what ensemble represents best thermal equilibrium in cores of globular 
clusters. However, that a stability limit exists is now commonly accepted. Globular clusters 



^■^Paper [g^ contains correct calculations but the wrong interpretation. The authors believed they were 
dealing with isolated systems when they dealt in fact with a grand canonical ensemble. 
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in quasi-thermal equilibrium evolve slowly towards greater entropy and density contrast due 
to stellar evaporation They then reach a limit beyond which they cannot stay isothermal. 
The great contribution of Antonov Q was to point out this instability. 

What happens when a core becomes unstable was explained by Lynden-Bell The 
explanation is now part of standard texts Stars diffuse towards the center which be- 
comes denser and denser while isothermal equilibrium is lost. This phenomena known as 
the gravothermal catastrophe leads eventually to core collapse. This has been confirmed by 
numerous studies mentioned in Meylan and Heggie |]50| . 

Observational evidence by various authors, see in particular Trager, Dorjovsky and King's 
paper |^5|, came in the late 1980's when CCD observations allowed a systematic investigation 
of the inner surface brightness profile of globular clusters. These authors classified globular 
clusters into two different classes. About 80% with a projected density profile that fitted 
Michie-King models and the rest with a density profile corresponding to a singular density 
in These cores are considered to have collapsed. 

The work of Antonov and Lynden-Bell is the single most important contribution to ob- 
servational astronomy based on thermodynamics of self-gravitating systems. 



(ii) Maclaurin and Jacobi ellipsoids. 



Liquid ellipsoids like isothermal spheres belong to standard literature. Chandrasekhar 
|To| devoted a whole book to the subject. Liquid ellipsoids have few independent variables 
and plenty of bifurcations. We give here one example in which a bifurcation is lifted with a 
small perturbation of the potential energy as discussed in Section 3. Solutions of Poisson's 
equation of equilibrium which involve elliptic integrals are taken from Chandrasekhar's book. 
The stability analysis presented here is taken from |34|. 



(a) Equilibrium configurations. 



Consider a self-gravitating ellipsoid with uniform density p and semi-axis a > b > c. It 
rotates around the c-axis with uniform and constant angular velocity T. The total mass M 
and angular momentum L are constant: 

M = —pabc , L = -Mfa^ + b'^)T = IT. (5.5) 
3 5 

The motion described in comoving coordinates along o, b, c has an effective potential energy 
V which is, see for instance |45] : 

^ = ^/pf^d'^ + |^- (5.6) 

U is the gravitational potential. There are three variables in V, a, b, c, and one constraint 
M. There are thus two independent variables, say, 

c c 

The explicit form of V is given in Chandrasekhar; in our notations, 

^ ' Jo y^il + iy){x^ + l^)ix^ + 1^) Xl+x2 ^ ^ 
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Equilibrium configurations X^{s) and X'^{s) are solutions of 

dV dV 

There exists one class of solutions to these equations for s < 0.769: the Maclaurin spheroids 
with a = b and an eccentricity e = (1 — < 0.813. For s — > 0, | = | ^ 1; spheroids 
turn into spheres of infinite radius and zero mass. When s > 0.769 there are two classes of 
solutions: Maclaurin spheroids with e > 0.813 and Jacobi ellipsoids with non equal axis. At 
s = 0.769 there is thus a bifurcation. 

The solutions are shown and discussed in details in Chandrasekhar's book together with 
their stability. Here we use the topology of linear series with conjugate variables with respect 
to minus the potential energy to find stability conditions. 

(b) Stability conditions. 

Let Fe{s) represent the extremal values of F. Fe{s) is equal to minus the potential energy 
times a constant. Equilibrium configurations are thus stable if Fe{s) is a maximum of F. 
The conjugate parameter of s with respect to Fe{s) is 



The linear series reproduced from |34| is shown in figure 4. We see the bifurcation appearing 
at point B where s = 0.769. Since homogeneous spheres are stable p6[ | Maclaurin spheroids 
are stable for < s < 0.769 and < e < 0.813. A change of stability can only occur at point 
B as we know but we cannot tell what happens at a bifurcation point just by looking at the 
linear series. 



Let us then look as suggested in (3.13) at solutions of modified equations like, for instance, 



dV dV 

with e small. The effect of e is to break the symmetry and lift the bifurcation. Figure 4 
represents also the linear series K{s) for the perturbed solutions with e = —0.01. Any other 
small value of e would have a similar effect on the topology but curves with e > cut those 
with e = and the drawing is not so nice. We now see a continuous line (sequence 1) 
which connects smoothly for e — > to the linear series of stable Maclaurin configurations 
for s < 0.769 and to Jacobi ellipsoids which exists only for s > 0.769. Sequence 1 has 
no vertical tangent. Thus by continuity Jacobi ellipsoids must be stable. Sequence 2 has 
a vertical tangent at point C. The CD branch identifies with the stable sequence 1 when 
e — > 0, and the CE branch represents necessarily unstable configurations with one unstable 
modef^ Taking e ^ 0, we conclude that Maclaurin spheroids with eccentricities e > 0.813 
are unstable. 

The stability limits have been known since Poincare's time [^]. The instability of Maclau- 
rin spheroids for e > 0.813 is secular; it shows up with a small viscous dissipation in a 
dynamical perturbation calculation [|l0| ]. 

(in) Final comments and more examples. 
^''The two Poincare coefficients of stability are different at point B, see 134]. 
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(a) Final Comments 



Thermodynamics of self-gravitational systems is helpful for calculating statistical equilib- 
rium configurations, thermodynamic functions, stability limits and the effect of fluctuations 
near instability. 

Globular clusters and other astronomical objects are only close to statistical equilibrium 
and evolve more often slowly. When evolution is slow enough these objects pass through 
a series of quasi-equilibrium configurations with ever increasing entropy. In this respect, 
linear series also mimic evolutionary tracks. Turning points indicate the limits beyond which 
evolution proceeds in a non isothermal way, generally towards higher central densities and 
core collapse. The fluctuations move somewhat the stability limits to lower temperatures. 

What happens next, once thermal equilibrium no more exists, is of immense interest but 
does not belong to this paper. Thermodynamics is a small chapter in theoretical astronomy 
though an interesting one. 



(b) More applications 



Thermodynamic methods as described in the present work have been used in various 
studies which we shall now mention briefly. One dimensional isothermal parallel sheets have 
been shown to be stable to one dimensional perturbations |^^. Isothermal axially-symmetric 
equilibria of gravitating rods are generally stable with respect to cylindrically symmetric 
perturbations ||3^. There is no gravothermal catastrophe as in isothermal spheres but no 
equilibria exist below some finite temperature. In contact with a heat bath slightly colder 
than that, the system collapses slowly, giving up an unlimited amount of energy. The effect 
of a mass spectrum on the stability limits of isothermal spheres has also been analyzed |38]. 

Isothermal spheres with a Fermi-Dirac distribution behave in some ways like systems 
of hard spheres particles: they have phase transitions from gaseous to core-halo structure. 
Stability limits with a Fermi-Dirac distribution were analyzed by Chavanis and Sommeria 
]ri| , see also [12|. A cutoff in potential has a similar effect Another example in which 
there may be phase transitions from white dwarfs to neutron stars but where the turning 
point method has not been used, though it would have proven extremely useful, is the study 
by Harrison, Thorne Wakano and Wheeler [21| of cold catalyzed matter. They minimize the 
mass-energy M. The only independent parameter is the number of baryons N . A plot of 
the conjugate parameter d{—M)/dN versus using their own table on page 152-153 would 
have shown the stability conditions and in particular the number of unstable modes in one 
stroke |37] as is clear from figure 5. 

Chavanis |14] also studied isothermal spheres in Newtonian mechanics but with a rela- 
tivistic equation of state. He found that instability sets in at smaller density contrasts but 
stronger binding energies. 

Stability conditions in relativistic spheres with energy cutoffs [2^26] have also not been 
treated by the turning point method but it is obvious from figures in |2S] that the their results 
might have been obtained in a more simple way. 

Since Bekenstein [|5|[P attributed an entropy to black holes and Hawking |22] found that 
they emit black body radiation of quantum origin, studies of their stability in various sur- 
roundings and with various charges have flourishes. After the seminal papers by Hawking 
|23| and by Gibbons and Perry|19] on the thermodynamics of black holes, there appeared 
a series of papers based on linear series devoted to Kerr black holes [E8| and Kerr-Newman 
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black holes [39|. A more detailed analysis of black holes in a cavity which takes account 



of the effect of fluctuations started in |53| was analyzed in greater details in ||5^. All those 
papers neglect the mass-energy of the cavity, usually a thin shell, that is supposed to play the 
role of a massless reflective wall. It has been pointed out that thin shells which do not 
allow infinite tensionsj^ must contain at least 30% of the total mass-energy. This may affect 
considerably the stability limits and throws some doubts on considerations, with sometimes 
far reaching consequences, that were derived from black holes in equilibrium with radiation 
in weightless cavities [|6l]] . 
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Appendices 



A Proof that ( p^.2c{| ) can be transformed into (^T2^) 



We start from (|2.23| ) with H = Hk + Hp defined in ( 2.21 ). The trick consists in introducing 
first a "continuous density" p(r): 



p{r) = ^ m6{r — rj). 



(A.l) 



The potential energy Hp can then be written 

Hp 



1 f Gp(r)P(r'),3,,3,, 



„ „ . „ . . (A.2) 
2Jv \r-r'\ ^ ' 

Next we discretize the whole space, cutting it into small cubes with an indice a and volume 



. Replacing p{r)d'^r in ( A.2 ) in each cube with indices a by p°'K^ the potential energy takes 
the form 

^P = ^EwV- (A.3) 



The matrix Uab is the discrete inverse of the Laplacian and more precisely 

A 



lim 

a,b-^ continuum 



AttG 



5{v-v'). 



(A.4) 



Let Wa represent the mean field in discretized form. Define a set of variables Xa to be 
integrated out soon: 

Xa = Wa- fi^'^Uabp''. (A.5) 

The following identityf^ is readily constructed with ([A.3| ) and ( |A.5|) by elementary algebraic 
manipulations: 



- pHp + hi-^'^^XaX,, = -p^/^Wap'' + ^U-^'^'WaWb. 



(A.6) 



Then taking the exponent of both sides of this expression and integrating over the whole 
domain of variation of the WaS and X^s leads to: 



-l3Hp 



[det(-27rua6)] 



-T / ^Y[dWa. (A.7) 

lTrii.„i,] \ 2 J— CO 



We now return to continuity, taking into account formula ( A.2 ) and ( A.4 ); this lets us write 
g-l3Hp ^Yie following form: 



vw, 



where 



B= lim J det{-2TrUab) and VW = lim TT'^^a- 

a,6^continuum ' a— ^continuum 



(A.8) 
(A.9) 



^This old identity due to Stratonovich ^] has proved useful in many-body physics. 
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B is certainly divergent; Uab tends to l/|r — r'|. We assume that a short distance cutoff 
will make B convergent. The nature of the cutoff is unimportant in our approximation - see 
section l(iii)(c). The expression for e~^^^ obtained in ( A.8 ) we put back into Q{E) in ( 2.23 ). 
It takes now little work to see that Q can be written as follows 



n 



1 



{2Tri)B Jb-i 



6+ioo /■+00 



-too J —oo 



where 



N 



is of the form ^ = — . 

NV 



(A.IO) 



(A.ll) 



This is the of the same form as ^ in (^^) . We may thus replace ^ with the same approximate 
expression as the one obtained for ^! in (|2.6|) to (p.9D : 



1 



in which 



N 



1 



1 



Ni a(N'-N) 



.N'=0 



da 



1 



a+ioo 



1 



^ N'\ 



Nf 



.N'=0 



5i-p2+m/3"^iy(r) 



We can also write ^ like this 



1 

27r« J a 



a+ioo 



(A.12) 
(A.13) 

(A.14) 



replacing \& by this expression in ( |A.10 ) gives n{E) as written in equation ( 2.2g| ) with a 
shown in equation ( p. 29 ). 

We followed closely the calculationsP^ in [p^ ] and |25| . We left aside all sorts of intricacies 
about existence and convergence problems in functional integrations. These were given more 
carefully considerations in |24| with appropriate references. 



B Calculation of fl{E) to order two in 6W - Formula ( \^.'6^ ) 

We start from a{a, f3, W) given in equation ( p. 29 ) which we expand in a Taylor series to order 
two in 5a = a — Oe, SP = (3 — j3e-,5W = W — We around their extremal values ae,Pe,We 
defined by ([23c|)-(|232|). There is no reason a priori to limit 6W to continuous perturbations 
except for the fact |2C] that non-continuous functions form a subset of measure zero. Here 
are the first order derivatives: 



da 
dW 



da 

da 

6Wd^r 



-N + / fduj 



da 
dp 



E- 



mp2 / 5Wfdco + 



1 



4ttG J- 



2m 2 



5WAWd-^r. 



fduj, 



(B.l) 



It will be of great help to define with an over-bar the mean values in the phase space of one 
particle. For instance: 



X 



1 

N 



XfeduJ 



(B.2) 



^''Paper [g^ contains correct calculations but the wrong interpretation. The authors believed they were 
dealing with isolated systems when they dealt in fact with a grand canonical ensemble. 
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Thus, the mean value of the total energy, equation ( p.32| ), can be written as follows: 

" 3 



E = Ek + Ep = N 



(PL] 



+ -mNU = N 



^ H — mil 
2/5e 2 



^ = 1 + lm(3^U. (B.3) 



\2mJ 2 

This expression will soon be useful. Another useful relation is the virial equality: 

N 2 [mN/V) Pmean 



(B.4) 



Here now are the second derivatives of a with respect to a, /?, W calculated at the point 
of extremum: 
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With these derivatives we may write the expression for 6'^ a at the extremum: 
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We now use ( |B.3D and (|B.4| ) to simplify some of the terms. Our aim is to reduce to a sum 
of squares that can be separately integrated in a, (5 and W or equivalently in 6a, 5(3 and 6W. 
We assume as usual that the exponent is steep and that order two is a good approximation. 
6a and SP vary along a line parallel to the imaginary axis between ±00. 
First we write 
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Second, 
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With these two results we can now see that the quadratic sum reduces to this: 
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The sum of 1' and 3', it is a difference of positive quantities: 
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We now observe that = eq.( |Z?.10| ) + eq.( B.ll ) + Terms(6 + 7). Let us introduce a couple 
of real variables which vary between ±00, similar to those introduced in equations ( |2.17| ): 
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In terms of a*, P* and 6W, 6"^ a reduces to 
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With this result we can calculate w to order two starting from ( |2.28| ), ( p.29| ); inserting 
£7 = 5 + 5^a into ( 2.28 ) we obtain 
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This is equivalent to formula ( 2.35 ). 
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It is worth noticing that 5^c7, see ( |B.9D , is equal to terms (6 + 7) minus two squares, see 
(^l0|) and pHI) . This imphes that: 



1 r+oo /A \ 

(5V < Terms (6 + 7) = - / 5W { —— + mpeP 5Wd?r = d'^acc (B.15) 

2 j-oo V47rG / 

is the second order term of the exponent for the grand canonical ensemble ace in- 
troduced in equation ( p.39D . The spectrum of eigenvalues of this quadratic form has been 
studied in detail in |2^] for isothermal spheres where it was shown that in a stable system 
^"^(^GC ~^ —oo. Thus 5^a — > — oo as well for such configurations. See section l(iii)(c) about 
this divergence. 



Figure captions 



Figure 1: This figure illustrates the point made in section 4 subsection (i) where full 
details are given. 

Figure 2: The figure represents the pair of conjugate thermodynamic functions P{E) in 
appropriate units for isothermal spheres described in section 4(i){a). The linear series gives 
stability limits in either isolated spheres or spheres in a heat bath. The parameter along the 
line, at vertical and at horizontal tangents is the density contrast TZ = p{center) / p[boundary) . 

Figure 3a and 3b: These figures represent the same pair of conjugate thermodynamic 
functions f3{E) as in figure 2 with different units (defined in the text) and is appropriate to 
detect stability limits of isothermal spheres in a grand canonical or a grand micro-canonical 
ensemble as described in section 4:(i){(3). Figure 3a is for —0.008 < < 0. Figure 3b is for 
< E < 0.6. For density contrasts 7^ ^ 1, — > and /3 oo. Thus the line with density 
contrasts TZ > 1.66 will, before E turns negative, cuts the line that comes from infinity . 

Figure 4: This figure represents the pair of conjugate variables for Maclaurin and Jacobi 
ellipsoids described in Afn). The figure is reproduced from |34] in which our small s was 
represented by S. S here represents the entropy. 

Figure 5: This is a topologically correct plot of d{—E)/ dN as a function of the number of 
baryons in units of Nq = 10^^. The parameter along the line is the central density in gr-cm~^ . 
The linear series shows six consecutive turning points calculated by Harrison. There are two 
stable branches corresponding respectively to cold white dwarfs and to neutron stars. More 
details are given in |21|. 
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